clc
clear

p = 100000;
rho = 0.1;
z = 80000;
a = 6371229;
omega = 7.292E-5;
gravity = 9.8;

nc_file = 'E:\Study\Models\FVM\run\output.nc';
u = ncread(nc_file,'u');
v = ncread(nc_file,'v');

dx = 1;

D2R = pi/180;
R2D = 180/pi;
dx = dx * D2R;
x_min = -pi/4+dx/2;
x_max =  pi/4-dx/2;
y_min = -pi/4+dx/2;
y_max =  pi/4-dx/2;
alpha = x_min:dx:x_max;
beta  = y_min:dx:y_max;

[alpha,beta] = meshgrid(alpha,beta);

X = tan(alpha);
Y = tan(beta);
sigma = sqrt(1+X.^2+Y.^2);

r = a + z;

SM3_1 = 2 * p / r;
SM3_2 = r * ( 1 + X.^2 ) .* ( 1 + Y.^2 ) ./ sigma.^4 .* ( ( 1 + X.^2 ) .* rho .* u .* u + ( 1 + Y.^2 ) .* rho .* v .* v - 2 .* X .* Y .* rho .* u .*v );
SM3 = SM3_1 + SM3_2;

FC3 = 2 * omega * r ./ sigma.^2 .* ( 1 + X.^2 ) .* rho .* u;

S3 = rho * gravity * a.^2 / r.^2;

S3 = S3 + SM3 + FC3;